library(tidyverse)
library(conformalInference.fd)
library(roahd)
library(lubridate)
library(mgcv)
library(mgcViz)
load("~/Documents/Nonparametric Statisics/Project/clean data/full_collisions.RData")
glimpse(full_collisions)
Rows: 2,585,717
Columns: 37
$ accident_index                              <chr> "200501BS00001", "200501BS00002", "200501BS0…
$ accident_year                               <dbl> 2005, 2005, 2005, 2005, 2005, 2005, 2005, 20…
$ accident_reference                          <chr> "01BS00001", "01BS00002", "01BS00003", "01BS…
$ location_easting_osgr                       <dbl> 525680, 524170, 524520, 526900, 528060, 5247…
$ location_northing_osgr                      <dbl> 178240, 181650, 182240, 177530, 179040, 1811…
$ longitude                                   <dbl> -0.191170, -0.211708, -0.206458, -0.173862, …
$ latitude                                    <dbl> 51.48910, 51.52007, 51.52530, 51.48244, 51.4…
$ police_force                                <fct> Metropolitan Police, Metropolitan Police, Me…
$ accident_severity                           <fct> Serious, Slight, Slight, Slight, Slight, Sli…
$ number_of_vehicles                          <dbl> 1, 1, 2, 1, 1, 2, 2, 1, 2, 2, 1, 2, 1, 2, 1,…
$ number_of_casualties                        <dbl> 1, 1, 1, 1, 1, 1, 1, 2, 2, 5, 1, 1, 1, 1, 2,…
$ date                                        <date> 2005-01-04, 2005-01-05, 2005-01-06, 2005-01…
$ day_of_week                                 <fct> Tuesday, Wednesday, Thursday, Friday, Monday…
$ time                                        <time> 17:42:00, 17:36:00, 00:15:00, 10:35:00, 21:…
$ local_authority_district                    <fct> "Kensington and Chelsea", "Kensington and Ch…
$ local_authority_ons_district                <chr> "E09000020", "E09000020", "E09000020", "E090…
$ local_authority_highway                     <chr> "E09000020", "E09000020", "E09000020", "E090…
$ first_road_class                            <fct> A, B, C, A, Unclassified, Unclassified, C, A…
$ first_road_number                           <dbl> 3218, 450, 0, 3220, 0, 0, 0, 315, 3212, 450,…
$ road_type                                   <fct> Single carriageway, Dual carriageway, Single…
$ speed_limit                                 <dbl> 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, …
$ junction_detail                             <fct> Not at junction or within 20 metres, Crossro…
$ junction_control                            <fct> Data missing or out of range, Auto traffic s…
$ second_road_class                           <fct> Not at junction or within 20 metres, C, Not …
$ second_road_number                          <dbl> -1, 0, -1, -1, -1, -1, 0, -1, 304, 0, 325, 3…
$ pedestrian_crossing_human_control           <fct> None within 50 metres, None within 50 metres…
$ pedestrian_crossing_physical_facilities     <fct> "Zebra", "Pedestrian phase at traffic signal…
$ light_conditions                            <fct> Daylight, Darkness - lights lit, Darkness - …
$ weather_conditions                          <fct> Raining no high winds, Fine no high winds, F…
$ road_surface_conditions                     <fct> Wet or damp, Dry, Dry, Dry, Wet or damp, Wet…
$ special_conditions_at_site                  <fct> None, None, None, None, None, Oil or diesel,…
$ carriageway_hazards                         <fct> None, None, None, None, None, None, None, No…
$ urban_or_rural_area                         <fct> Urban, Urban, Urban, Urban, Urban, Urban, Ur…
$ did_police_officer_attend_scene_of_accident <fct> Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes,…
$ trunk_road_flag                             <fct> Non-trunk, Non-trunk, Non-trunk, Non-trunk, …
$ lsoa_of_accident_location                   <chr> "E01002849", "E01002909", "E01002857", "E010…
$ datetime                                    <dttm> 2005-01-04 17:42:00, 2005-01-05 17:36:00, 2…

we need to create the dataset in a way that they cotain all of the days and a zero if no crashes happened in that day:

all_dates <- unique(full_collisions$date)
all_police_forces <- unique(full_collisions$police_force)

full_combinations <- expand.grid(date = all_dates, police_force = all_police_forces)

crashes_per_day_police <- full_combinations %>%
  left_join(full_collisions %>% group_by(date, police_force) %>% summarize(number_of_crashes = n()),
            by = c("date", "police_force")) %>%
  replace_na(list(number_of_crashes = 0))
`summarise()` has grouped output by 'date'. You can override using the `.groups` argument.
glimpse(crashes_per_day_police)
Rows: 341,848
Columns: 3
$ date              <date> 2005-01-04, 2005-01-05, 2005-01-06, 2005-01-07, 2005-01-10, 2005-01-1…
$ police_force      <fct> Metropolitan Police, Metropolitan Police, Metropolitan Police, Metropo…
$ number_of_crashes <int> 58, 63, 64, 59, 67, 76, 86, 82, 57, 57, 85, 92, 82, 79, 64, 80, 78, 80…
crashes_per_day_police <- crashes_per_day_police %>% mutate(day_of_year = yday(date),
                                                            day = factor(wday(date)),year = year(date))

removing the 29 of february for dimensioality problems in the vectors below

crashes_per_day_police <- crashes_per_day_police %>% filter(!(month(date)==2 & day(date)==29))
train_years <- ceiling(seq(2005,2021,by = 1.5))
cal_years <- seq(2006,2021,3)

df_train <- crashes_per_day_police %>% filter(year %in% train_years) 
df_cal <- crashes_per_day_police %>% filter(year %in% cal_years) 
df_test <- crashes_per_day_police %>% filter(year == 2022)

starting with a gam and mixed effects

gam_train <- bam(number_of_crashes ~ day + year +
                                  s(day_of_year,k = 53, bs = "cr") + 
                       s(police_force, bs = "re"), 
                data = df_train, family=quasipoisson(), method='REML')
summary(gam_train)

Family: quasipoisson 
Link function: log 

Formula:
number_of_crashes ~ day + year + s(day_of_year, k = 53, bs = "cr") + 
    s(police_force, bs = "re")

Parametric coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 81.2959321  0.4499718  180.67   <2e-16 ***
day2         0.2543089  0.0040811   62.31   <2e-16 ***
day3         0.3038483  0.0040350   75.30   <2e-16 ***
day4         0.3109257  0.0040289   77.17   <2e-16 ***
day5         0.3146008  0.0040276   78.11   <2e-16 ***
day6         0.3942311  0.0039628   99.48   <2e-16 ***
day7         0.1876563  0.0041400   45.33   <2e-16 ***
year        -0.0396915  0.0002158 -183.97   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                  edf Ref.df       F p-value    
s(day_of_year)  49.18   51.6   157.6  <2e-16 ***
s(police_force) 50.99   51.0 16496.4  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.851   Deviance explained = 80.6%
-REML = 3.5021e+05  Scale est. = 1.6704    n = 208780
b <- getViz(gam_train)
print(plot(b, allTerms = T), pages = 1)

1 is sunday here

newd <- with(df_train, data.frame(day_of_year = 120,
                                  day = 1,
                                  year = 2022,
                                  police_force = levels(police_force)))

p <- predict(gam_train, newd, type = "terms", se.fit = TRUE)
re <- p[["fit"]][ , "s(police_force)"]
se <- p[["se.fit"]][ , "s(police_force)"]

data <- data.frame(police_force = levels(df_train$police_force), effect = re,std.err = se)
data %>% 
ggplot(aes(effect, fct_reorder(police_force, effect))) +
  geom_vline(xintercept = 0, color = "gray50", lty = 2, linewidth = 1.2) +
  geom_errorbar(aes(
    xmin = effect - 1.96*std.err,
    xmax = effect + 1.96*std.err),width = 0.5, alpha = 0.7) +
  geom_point(size = 1,colour = "blue") + theme_minimal() + 
  labs(y = "District", x = "Random effect") + ggtitle("Random effects for each police force") + 
  theme(legend.position = "none",plot.title = element_text(size = 16,hjust = 0.5)) 

check(b)

we now need to compute the bands in a functional case

res_train <- cbind(df_train,preds = predict(gam_train, df_train, type = "response")) %>%
  arrange(police_force,date) %>% mutate(absdiff = abs(number_of_crashes-preds))

res_cal <- cbind(df_cal,preds = predict(gam_train, df_cal, type = "response")) %>%
  arrange(police_force,date) %>% mutate(absdiff = abs(number_of_crashes-preds))

res_test <- cbind(df_test,preds = predict(gam_train, df_test, type = "response")) %>%
  arrange(police_force,date) %>% mutate(absdiff = abs(number_of_crashes-preds))
n_pol <- length(levels(crashes_per_day_police$police_force))

m = 12
l = 6
alpha = 0.1
S = matrix(nrow = n_pol,ncol = 365)
k = rep(0,n_pol)
x = 1:365
line_integral = function(x, y) {
  dx = diff(x)
  end = length(y)
  my = (y[1:(end - 1)] + y[2:end]) / 2
  sum(dx *my)
} 

ceiling((m+1)*(1-alpha))

we use the largest value for gamma H1 = I1

computing the modulation functions:

for(i in 1:n_pol){
  ad_train <- res_train %>% filter(police_force == levels(police_force)[i]) %>%
  dplyr::select(absdiff)
  a_mat <- matrix(ad_train$absdiff,byrow = F,nrow = 365)
  ncm = apply(a_mat,1,max)
  den = line_integral(x, ncm)
  s = ncm/den
  S[i,] = s
}

computing the k:

ceiling((l + 1)*(1-alpha))

we use the biggest value for the choice of k

for(i in 1:n_pol){
  ad_cal <- res_cal %>% filter(police_force == levels(police_force)[i]) %>%
  dplyr::select(absdiff)
  a_mat <- matrix(ad_cal$absdiff,byrow = F,nrow = 365)
  ncs = apply(a_mat/S[i,],1,max)
  k[i] = max(ncs)
}
i <- 50

df_filtered <- res_test %>% filter(police_force == levels(police_force)[i]) %>%
  arrange(day_of_year)

df_plot <- data.frame(n = df_filtered$number_of_crashes,pred = df_filtered$preds, 
                      ymax = df_filtered$preds + k[i]*S[i,],
                      ymin = pmax( df_filtered$preds - k[i]*S[i,], rep(0, 365)),td = 1:365)

df_plot %>% ggplot() + geom_ribbon(aes(x=td,y = pred,ymin = ymin, ymax = ymax), 
                                   fill = "blue", alpha = 0.25) + 
  geom_line(aes(x=td,y = n),linewidth=1,colour = "black") + 
  geom_line(aes(x=td,y = pred),linewidth=1,colour = "blue") + labs(x = "Day", y = "Number of accients",
                                                     title = levels(df_test$police_force)[i])

this is not great:

i <- 6

df_filtered <- res_test %>% filter(police_force == levels(police_force)[i]) %>%
  arrange(day_of_year)

df_plot <- data.frame(n = df_filtered$number_of_crashes,pred = df_filtered$preds, 
                      ymax = df_filtered$preds + k[i]*S[i,],
                      ymin = pmax( df_filtered$preds - k[i]*S[i,], rep(0, 365)),td = 1:365)

df_plot %>% ggplot() + geom_ribbon(aes(x=td,y = pred,ymin = ymin, ymax = ymax), 
                                   fill = "blue", alpha = 0.25) + 
  geom_line(aes(x=td,y = n),linewidth=1,colour = "black") + 
  geom_line(aes(x=td,y = pred),linewidth=1,colour = "blue") + labs(x = "Day", y = "Number of accients",
                                                     title = levels(df_test$police_force)[i])

i <- 30

df_filtered <- res_test %>% filter(police_force == levels(police_force)[i]) %>%
  arrange(day_of_year)

df_plot <- data.frame(n = df_filtered$number_of_crashes,pred = df_filtered$preds, 
                      ymax = df_filtered$preds + k[i]*S[i,],
                      ymin = pmax( df_filtered$preds - k[i]*S[i,], rep(0, 365)),td = 1:365)

df_plot %>% ggplot() + geom_ribbon(aes(x=td,y = pred,ymin = ymin, ymax = ymax), 
                                   fill = "blue", alpha = 0.25) + 
  geom_line(aes(x=td,y = n),linewidth=1,colour = "black") + 
  geom_line(aes(x=td,y = pred),linewidth=1,colour = "blue") + labs(x = "Day", y = "Number of accients",
                                                     title = levels(df_test$police_force)[i])

this are too big, it is due to the fact that we are using a mixd effects model, so metropolitan police will contaminate all of the data.

we can fix this by making a gam for each police force:

gams <- list()

S = matrix(nrow = n_pol,ncol = 365)
k = rep(0,n_pol)

for(i in 1:n_pol){
  
  df_train_mid <- df_train %>% filter(police_force == levels(police_force)[i])
  df_cal_mid <- df_cal %>% filter(police_force == levels(police_force)[i])
  
  gam <- bam(number_of_crashes ~ day + year +
                                  s(day_of_year,k = 53, bs = "cr"), 
                data = df_train_mid, family=quasipoisson(), method='REML')
  
  res_train_mid <- cbind(df_train_mid,preds = predict(gam, df_train_mid, type = "response")) %>%
  arrange(date) %>% mutate(absdiff = abs(number_of_crashes-preds))

  res_cal_mid <- cbind(df_cal_mid,preds = predict(gam, df_cal_mid, type = "response")) %>%
  arrange(date) %>% mutate(absdiff = abs(number_of_crashes-preds))
  
  
  ad_train <- res_train_mid %>% filter(police_force == levels(police_force)[i]) %>%
  dplyr::select(absdiff)
  a_mat <- matrix(ad_train$absdiff,byrow = F,nrow = 365)
  ncm = apply(a_mat,1,max)
  den = line_integral(x, ncm)
  s = ncm/den
  S[i,] = s
  ad_cal <- res_cal_mid %>% filter(police_force == levels(police_force)[i]) %>%
  dplyr::select(absdiff)
  a_mat <- matrix(ad_cal$absdiff,byrow = F,nrow = 365)
  ncs = apply(a_mat/S[i,],1,max)
  k[i] = max(ncs)
  
  gams[[i]] <- gam
}

midlands

i <- 50
df_test_filtered <- df_test %>% filter(police_force == levels(police_force)[i])
preds <- predict(gams[[i]], df_test_filtered, type = "response")
df_test_filtered  <- cbind(df_test_filtered,preds) %>% arrange(day_of_year)
df_plot <- data.frame(n = df_test_filtered$number_of_crashes,pred = df_test_filtered$preds, 
                      ymax = df_test_filtered$preds + k[i]*S[i,],
                      ymin = pmax( df_test_filtered$preds - k[i]*S[i,], rep(0, 365)),td = 1:365)

df_plot %>% ggplot() + geom_ribbon(aes(x=td,y = pred,ymin = ymin, ymax = ymax), 
                                   fill = "blue", alpha = 0.25) + 
  geom_line(aes(x=td,y = n),linewidth=1,colour = "black") + 
  geom_line(aes(x=td,y = pred),linewidth=1,colour = "blue") + labs(x = "Day", y = "Number of accients",
                                                     title = levels(df_test$police_force)[i])

London:

i <- 30
df_test_filtered <- df_test %>% filter(police_force == levels(police_force)[i])
preds <- predict(gams[[i]], df_test_filtered, type = "response")
df_test_filtered  <- cbind(df_test_filtered,preds) %>% arrange(day_of_year)
df_plot <- data.frame(n = df_test_filtered$number_of_crashes,pred = df_test_filtered$preds, 
                      ymax = df_test_filtered$preds + k[i]*S[i,],
                      ymin = pmax( df_test_filtered$preds - k[i]*S[i,], rep(0, 365)),td = 1:365)

df_plot %>% ggplot() + geom_ribbon(aes(x=td,y = pred,ymin = ymin, ymax = ymax), 
                                   fill = "blue", alpha = 0.25) + 
  geom_line(aes(x=td,y = n),linewidth=1,colour = "black") + 
  geom_line(aes(x=td,y = pred),linewidth=1,colour = "blue") + labs(x = "Day", y = "Number of accients",
                                                     title = levels(df_test$police_force)[i])

the model is better but still no great

we probably need more data.

using all the available data:

full_collisioin_data <- readr::read_csv(file="dft-road-casualty-statistics-collision-1979-latest-published-year.csv") 

spec(full_collisioin_data)
glimpse(full_collisioin_data)

full_collisions_all <- stats19::format_collisions(full_collisioin_data)

glimpse(full_collisions_all)
all_dates <- unique(full_collisions_all$date)
all_police_forces <- unique(full_collisions_all$police_force)

full_combinations <- expand.grid(date = all_dates, police_force = all_police_forces)

crashes_per_day_police <- full_combinations %>%
  left_join(full_collisions_all %>% group_by(date, police_force) %>% 
              summarize(number_of_crashes = n()),
            by = c("date", "police_force")) %>%
  replace_na(list(number_of_crashes = 0))
`summarise()` has grouped output by 'date'. You can override using the `.groups` argument.
glimpse(crashes_per_day_police)
Rows: 835,744
Columns: 3
$ date              <date> 1979-01-18, 1979-01-01, 1979-01-02, 1979-01-03, 1979-01-04, 1979-01-0…
$ police_force      <fct> Metropolitan Police, Metropolitan Police, Metropolitan Police, Metropo…
$ number_of_crashes <int> 115, 39, 51, 59, 55, 84, 96, 72, 114, 132, 168, 128, 211, 145, 103, 14…
crashes_per_day_police <- crashes_per_day_police %>% mutate(day_of_year = yday(date),
                                                            day = factor(wday(date)),year = year(date))

removing the 29 of february for dimensioality problems in the vectors below

crashes_per_day_police <- crashes_per_day_police %>% filter(!(month(date)==2 & day(date)==29))
train_years <- ceiling(seq(1979,2021,by = 1.5))
cal_years <- seq(1979,2021,3)

df_train <- crashes_per_day_police %>% filter(year %in% train_years) 
df_cal <- crashes_per_day_police %>% filter(year %in% cal_years) 
df_test <- crashes_per_day_police %>% filter(year == 2022)

starting with a gam and mixed effects

gam_train <- bam(number_of_crashes ~ day + year +
                                  s(day_of_year,k = 53, bs = "cr") + 
                       s(police_force, bs = "re"), 
                data = df_train, family=quasipoisson(), method='REML')
summary(gam_train)

Family: quasipoisson 
Link function: log 

Formula:
number_of_crashes ~ day + year + s(day_of_year, k = 53, bs = "cr") + 
    s(police_force, bs = "re")

Parametric coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.030e+01  1.526e-01   264.1   <2e-16 ***
day2         2.351e-01  2.333e-03   100.7   <2e-16 ***
day3         2.505e-01  2.326e-03   107.7   <2e-16 ***
day4         2.651e-01  2.319e-03   114.3   <2e-16 ***
day5         2.991e-01  2.301e-03   129.9   <2e-16 ***
day6         4.129e-01  2.248e-03   183.7   <2e-16 ***
day7         2.345e-01  2.334e-03   100.5   <2e-16 ***
year        -1.929e-02  4.717e-05  -409.1   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                  edf Ref.df       F p-value    
s(day_of_year)  51.09  51.96   508.8  <2e-16 ***
s(police_force) 51.00  51.00 54521.2  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.894   Deviance explained = 84.1%
-REML = 9.6883e+05  Scale est. = 1.9754    n = 550420
b <- getViz(gam_train)
print(plot(b, allTerms = T), pages = 1)

1 is sunday here

newd <- with(df_train, data.frame(day_of_year = 120,
                                  day = 1,
                                  year = 2022,
                                  police_force = levels(police_force)))

p <- predict(gam_train, newd, type = "terms", se.fit = TRUE)
re <- p[["fit"]][ , "s(police_force)"]
se <- p[["se.fit"]][ , "s(police_force)"]

data <- data.frame(police_force = levels(df_train$police_force), effect = re,std.err = se)
data %>% 
ggplot(aes(effect, fct_reorder(police_force, effect))) +
  geom_vline(xintercept = 0, color = "gray50", lty = 2, linewidth = 1.2) +
  geom_errorbar(aes(
    xmin = effect - 1.96*std.err,
    xmax = effect + 1.96*std.err),width = 0.5, alpha = 0.7) +
  geom_point(size = 1,colour = "blue") + theme_minimal() + 
  labs(y = "District", x = "Random effect") + ggtitle("Random effects for each police force") + 
  theme(legend.position = "none",plot.title = element_text(size = 16,hjust = 0.5)) 

check(b)

Method: REML   Optimizer: outer newton
full convergence after 7 iterations.
Gradient range [-0.1344576,0.130992]
(score 968827.9 & scale 1.975359).
Hessian positive definite, eigenvalue range [23.46106,275205.6].
Model rank =  112 / 112 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

                  k'  edf k-index p-value
s(day_of_year)  52.0 51.1    0.99    0.56
s(police_force) 52.0 51.0      NA      NA

we now need to compute the bands in a functional case

res_train <- cbind(df_train,preds = predict(gam_train, df_train, type = "response")) %>%
  arrange(police_force,date) %>% mutate(absdiff = abs(number_of_crashes-preds))

res_cal <- cbind(df_cal,preds = predict(gam_train, df_cal, type = "response")) %>%
  arrange(police_force,date) %>% mutate(absdiff = abs(number_of_crashes-preds))

res_test <- cbind(df_test,preds = predict(gam_train, df_test, type = "response")) %>%
  arrange(police_force,date) %>% mutate(absdiff = abs(number_of_crashes-preds))
n_pol <- length(levels(crashes_per_day_police$police_force))

m = length(train_years)
l = length(cal_years)
alpha = 0.1
S = matrix(nrow = n_pol,ncol = 365)
k = rep(0,n_pol)
x = 1:365
line_integral = function(x, y) {
  dx = diff(x)
  end = length(y)
  my = (y[1:(end - 1)] + y[2:end]) / 2
  sum(dx *my)
} 

ceiling((m+1)*(1-alpha))
[1] 27

H1 != I1 now

computing the modulation functions:

for(i in 1:n_pol){
  ad_train <- res_train %>% filter(police_force == levels(police_force)[i]) %>%
  dplyr::select(absdiff)
  a_mat <- matrix(ad_train$absdiff,byrow = F,nrow = 365)
  setcfr <- apply(a_mat,2,max)
  gamma <- sort(setcfr)[ceiling((m+1)*(1-alpha))]
  a_mat_filt <- a_mat[,setcfr<=gamma]
  ncm = apply(a_mat_filt,1,max)
  ncm_sort <- sort(ncm)
  den = line_integral(x, ncm)
  s = ncm/den
  S[i,] = s
}

computing the k:

ceiling((l + 1)*(1-alpha))
[1] 15

we use the biggest value for the choice of k

for(i in 1:n_pol){
  ad_cal <- res_cal %>% filter(police_force == levels(police_force)[i]) %>%
  dplyr::select(absdiff)
  a_mat <- matrix(ad_cal$absdiff,byrow = F,nrow = 365)
  ncs = apply(a_mat/S[i,],1,max)
  k[i] = max(ncs)
}
i <- 50

df_filtered <- res_test %>% filter(police_force == levels(police_force)[i]) %>%
  arrange(day_of_year)

df_plot <- data.frame(n = df_filtered$number_of_crashes,pred = df_filtered$preds, 
                      ymax = df_filtered$preds + k[i]*S[i,],
                      ymin = pmax( df_filtered$preds - k[i]*S[i,], rep(0, 365)),td = 1:365)

df_plot %>% ggplot() + geom_ribbon(aes(x=td,y = pred,ymin = ymin, ymax = ymax), 
                                   fill = "blue", alpha = 0.25) + 
  geom_line(aes(x=td,y = n),linewidth=1,colour = "black") + 
  geom_line(aes(x=td,y = pred),linewidth=1,colour = "blue") + labs(x = "Day", y = "Number of accients",
                                                     title = levels(df_test$police_force)[i])

this is not great:

i <- 6

df_filtered <- res_test %>% filter(police_force == levels(police_force)[i]) %>%
  arrange(day_of_year)

df_plot <- data.frame(n = df_filtered$number_of_crashes,pred = df_filtered$preds, 
                      ymax = df_filtered$preds + k[i]*S[i,],
                      ymin = pmax( df_filtered$preds - k[i]*S[i,], rep(0, 365)),td = 1:365)

df_plot %>% ggplot() + geom_ribbon(aes(x=td,y = pred,ymin = ymin, ymax = ymax), 
                                   fill = "blue", alpha = 0.25) + 
  geom_line(aes(x=td,y = n),linewidth=1,colour = "black") + 
  geom_line(aes(x=td,y = pred),linewidth=1,colour = "blue") + labs(x = "Day", y = "Number of accients",
                                                     title = levels(df_test$police_force)[i])

i <- 30

df_filtered <- res_test %>% filter(police_force == levels(police_force)[i]) %>%
  arrange(day_of_year)

df_plot <- data.frame(n = df_filtered$number_of_crashes,pred = df_filtered$preds, 
                      ymax = df_filtered$preds + k[i]*S[i,],
                      ymin = pmax( df_filtered$preds - k[i]*S[i,], rep(0, 365)),td = 1:365)

df_plot %>% ggplot() + geom_ribbon(aes(x=td,y = pred,ymin = ymin, ymax = ymax), 
                                   fill = "blue", alpha = 0.25) + 
  geom_line(aes(x=td,y = n),linewidth=1,colour = "black") + 
  geom_line(aes(x=td,y = pred),linewidth=1,colour = "blue") + labs(x = "Day", y = "Number of accients",
                                                     title = levels(df_test$police_force)[i])

this is eve worse tha before:

doing it for a single police force at a time:

gams <- list()

S = matrix(nrow = n_pol,ncol = 365)
k = rep(0,n_pol)

for(i in 1:n_pol){
  
  df_train_mid <- df_train %>% filter(police_force == levels(police_force)[i])
  df_cal_mid <- df_cal %>% filter(police_force == levels(police_force)[i])
  
  gam <- bam(number_of_crashes ~ day + year +
                                  s(day_of_year,k = 53, bs = "cr"), 
                data = df_train_mid, family=quasipoisson(), method='REML')
  
  res_train_mid <- cbind(df_train_mid,preds = predict(gam, df_train_mid, type = "response")) %>%
  arrange(date) %>% mutate(absdiff = abs(number_of_crashes-preds))

  res_cal_mid <- cbind(df_cal_mid,preds = predict(gam, df_cal_mid, type = "response")) %>%
  arrange(date) %>% mutate(absdiff = abs(number_of_crashes-preds))
  
  
  
  ad_train <- res_train_mid %>% filter(police_force == levels(police_force)[i]) %>%
  dplyr::select(absdiff)
  a_mat <- matrix(ad_train$absdiff,byrow = F,nrow = 365)
  setcfr <- apply(a_mat,2,max)
  gamma <- sort(setcfr)[ceiling((m+1)*(1-alpha))]
  a_mat_filt <- a_mat[,setcfr<=gamma]
  ncm = apply(a_mat_filt,1,max)
  ncm_sort <- sort(ncm)
  den = line_integral(x, ncm)
  s = ncm/den
  S[i,] = s
  
  ad_cal <- res_cal_mid %>% filter(police_force == levels(police_force)[i]) %>%
  dplyr::select(absdiff)
  a_mat <- matrix(ad_cal$absdiff,byrow = F,nrow = 365)
  ncs = apply(a_mat/S[i,],1,max)
  k[i] = max(ncs)
  
  gams[[i]] <- gam
}

midlands

i <- 50
df_test_filtered <- df_test %>% filter(police_force == levels(police_force)[i])
preds <- predict(gams[[i]], df_test_filtered, type = "response")
df_test_filtered  <- cbind(df_test_filtered,preds) %>% arrange(day_of_year)
df_plot <- data.frame(n = df_test_filtered$number_of_crashes,pred = df_test_filtered$preds, 
                      ymax = df_test_filtered$preds + k[i]*S[i,],
                      ymin = pmax( df_test_filtered$preds - k[i]*S[i,], rep(0, 365)),td = 1:365)

df_plot %>% ggplot() + geom_ribbon(aes(x=td,y = pred,ymin = ymin, ymax = ymax), 
                                   fill = "blue", alpha = 0.25) + 
  geom_line(aes(x=td,y = n),linewidth=1,colour = "black") + 
  geom_line(aes(x=td,y = pred),linewidth=1,colour = "blue") + labs(x = "Day", y = "Number of accients",
                                                     title = levels(df_test$police_force)[i])

London:

i <- 30
df_test_filtered <- df_test %>% filter(police_force == levels(police_force)[i])
preds <- predict(gams[[i]], df_test_filtered, type = "response")
df_test_filtered  <- cbind(df_test_filtered,preds) %>% arrange(day_of_year)
df_plot <- data.frame(n = df_test_filtered$number_of_crashes,pred = df_test_filtered$preds, 
                      ymax = df_test_filtered$preds + k[i]*S[i,],
                      ymin = pmax( df_test_filtered$preds - k[i]*S[i,], rep(0, 365)),td = 1:365)

df_plot %>% ggplot() + geom_ribbon(aes(x=td,y = pred,ymin = ymin, ymax = ymax), 
                                   fill = "blue", alpha = 0.25) + 
  geom_line(aes(x=td,y = n),linewidth=1,colour = "black") + 
  geom_line(aes(x=td,y = pred),linewidth=1,colour = "blue") + labs(x = "Day", y = "Number of accients",
                                                     title = levels(df_test$police_force)[i])

the model is still ot graet

LS0tCnRpdGxlOiAiTm9ucGFyYW1ldHJpYyBBbmFseXNpcyBvZiBVSyBSb2FkIEFjY2lkZW50cyIKc3VidGl0bGU6ICJmdW5jdGlvbmFsIGNwIGZvciBlYWNoIHBvbGljZSBmb3JjZSIKYXV0aG9yOgogICAgLSAiVmFsZXJpYSBJYXBhb2xvIgogICAgLSAiT3N3YWxkbyBNb3JhbGVzIgogICAgLSAiUmljY2FyZG8gTW9yYW5kaSIKICAgIC0gIkFieWxhaSBPcnluYmFzc2FyIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpgYGB7ciBzZXR1cCwgZWNobyA9IEZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQoCiAgICBlY2hvID0gVFJVRSwKICAgICNkZXYgPSBjKCdwZGYnKSwKICAgIGZpZy5hbGlnbiA9ICdjZW50ZXInLAogICAgI2ZpZy5wYXRoID0gJ291dHB1dC8nLAogICAgZmlnLmhlaWdodCA9IDYsCiAgICBmaWcud2lkdGggPSAxMgopCmBgYAoKYGBge3J9CmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KGNvbmZvcm1hbEluZmVyZW5jZS5mZCkKbGlicmFyeShyb2FoZCkKbGlicmFyeShsdWJyaWRhdGUpCmxpYnJhcnkobWdjdikKbGlicmFyeShtZ2NWaXopCmBgYAoKYGBge3J9CmxvYWQoIn4vRG9jdW1lbnRzL05vbnBhcmFtZXRyaWMgU3RhdGlzaWNzL1Byb2plY3QvY2xlYW4gZGF0YS9mdWxsX2NvbGxpc2lvbnMuUkRhdGEiKQpgYGAKCmBgYHtyfQpnbGltcHNlKGZ1bGxfY29sbGlzaW9ucykKYGBgCgp3ZSBuZWVkIHRvIGNyZWF0ZSB0aGUgZGF0YXNldCBpbiBhIHdheSB0aGF0IHRoZXkgY290YWluIGFsbCBvZiB0aGUgZGF5cyBhbmQgYSB6ZXJvIGlmIG5vIGNyYXNoZXMgaGFwcGVuZWQgaW4gdGhhdCBkYXk6CgpgYGB7cn0KYWxsX2RhdGVzIDwtIHVuaXF1ZShmdWxsX2NvbGxpc2lvbnMkZGF0ZSkKYWxsX3BvbGljZV9mb3JjZXMgPC0gdW5pcXVlKGZ1bGxfY29sbGlzaW9ucyRwb2xpY2VfZm9yY2UpCgpmdWxsX2NvbWJpbmF0aW9ucyA8LSBleHBhbmQuZ3JpZChkYXRlID0gYWxsX2RhdGVzLCBwb2xpY2VfZm9yY2UgPSBhbGxfcG9saWNlX2ZvcmNlcykKCmNyYXNoZXNfcGVyX2RheV9wb2xpY2UgPC0gZnVsbF9jb21iaW5hdGlvbnMgJT4lCiAgbGVmdF9qb2luKGZ1bGxfY29sbGlzaW9ucyAlPiUgZ3JvdXBfYnkoZGF0ZSwgcG9saWNlX2ZvcmNlKSAlPiUgc3VtbWFyaXplKG51bWJlcl9vZl9jcmFzaGVzID0gbigpKSwKICAgICAgICAgICAgYnkgPSBjKCJkYXRlIiwgInBvbGljZV9mb3JjZSIpKSAlPiUKICByZXBsYWNlX25hKGxpc3QobnVtYmVyX29mX2NyYXNoZXMgPSAwKSkKCmdsaW1wc2UoY3Jhc2hlc19wZXJfZGF5X3BvbGljZSkKYGBgCgpgYGB7cn0KY3Jhc2hlc19wZXJfZGF5X3BvbGljZSA8LSBjcmFzaGVzX3Blcl9kYXlfcG9saWNlICU+JSBtdXRhdGUoZGF5X29mX3llYXIgPSB5ZGF5KGRhdGUpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBkYXkgPSBmYWN0b3Iod2RheShkYXRlKSkseWVhciA9IHllYXIoZGF0ZSkpCmBgYAoKcmVtb3ZpbmcgdGhlIDI5IG9mIGZlYnJ1YXJ5IGZvciBkaW1lbnNpb2FsaXR5IHByb2JsZW1zIGluIHRoZSB2ZWN0b3JzIGJlbG93CgpgYGB7cn0KY3Jhc2hlc19wZXJfZGF5X3BvbGljZSA8LSBjcmFzaGVzX3Blcl9kYXlfcG9saWNlICU+JSBmaWx0ZXIoIShtb250aChkYXRlKT09MiAmIGRheShkYXRlKT09MjkpKQpgYGAKCmBgYHtyfQp0cmFpbl95ZWFycyA8LSBjZWlsaW5nKHNlcSgyMDA1LDIwMjEsYnkgPSAxLjUpKQpjYWxfeWVhcnMgPC0gc2VxKDIwMDYsMjAyMSwzKQoKZGZfdHJhaW4gPC0gY3Jhc2hlc19wZXJfZGF5X3BvbGljZSAlPiUgZmlsdGVyKHllYXIgJWluJSB0cmFpbl95ZWFycykgCmRmX2NhbCA8LSBjcmFzaGVzX3Blcl9kYXlfcG9saWNlICU+JSBmaWx0ZXIoeWVhciAlaW4lIGNhbF95ZWFycykgCmRmX3Rlc3QgPC0gY3Jhc2hlc19wZXJfZGF5X3BvbGljZSAlPiUgZmlsdGVyKHllYXIgPT0gMjAyMikKYGBgCgpzdGFydGluZyB3aXRoIGEgZ2FtIGFuZCBtaXhlZCBlZmZlY3RzCgpgYGB7cn0KZ2FtX3RyYWluIDwtIGJhbShudW1iZXJfb2ZfY3Jhc2hlcyB+IGRheSArIHllYXIgKwogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcyhkYXlfb2ZfeWVhcixrID0gNTMsIGJzID0gImNyIikgKyAKICAgICAgICAgICAgICAgICAgICAgICBzKHBvbGljZV9mb3JjZSwgYnMgPSAicmUiKSwgCiAgICAgICAgICAgICAgICBkYXRhID0gZGZfdHJhaW4sIGZhbWlseT1xdWFzaXBvaXNzb24oKSwgbWV0aG9kPSdSRU1MJykKYGBgCgpgYGB7cn0Kc3VtbWFyeShnYW1fdHJhaW4pCmBgYAoKYGBge3J9CmIgPC0gZ2V0Vml6KGdhbV90cmFpbikKYGBgCgpgYGB7cn0KcHJpbnQocGxvdChiLCBhbGxUZXJtcyA9IFQpLCBwYWdlcyA9IDEpCmBgYAoKMSBpcyBzdW5kYXkgaGVyZQoKYGBge3J9Cm5ld2QgPC0gd2l0aChkZl90cmFpbiwgZGF0YS5mcmFtZShkYXlfb2ZfeWVhciA9IDEyMCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGRheSA9IDEsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5ZWFyID0gMjAyMiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHBvbGljZV9mb3JjZSA9IGxldmVscyhwb2xpY2VfZm9yY2UpKSkKCnAgPC0gcHJlZGljdChnYW1fdHJhaW4sIG5ld2QsIHR5cGUgPSAidGVybXMiLCBzZS5maXQgPSBUUlVFKQpyZSA8LSBwW1siZml0Il1dWyAsICJzKHBvbGljZV9mb3JjZSkiXQpzZSA8LSBwW1sic2UuZml0Il1dWyAsICJzKHBvbGljZV9mb3JjZSkiXQoKZGF0YSA8LSBkYXRhLmZyYW1lKHBvbGljZV9mb3JjZSA9IGxldmVscyhkZl90cmFpbiRwb2xpY2VfZm9yY2UpLCBlZmZlY3QgPSByZSxzdGQuZXJyID0gc2UpCmBgYAoKYGBge3J9CmRhdGEgJT4lIApnZ3Bsb3QoYWVzKGVmZmVjdCwgZmN0X3Jlb3JkZXIocG9saWNlX2ZvcmNlLCBlZmZlY3QpKSkgKwogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IDAsIGNvbG9yID0gImdyYXk1MCIsIGx0eSA9IDIsIGxpbmV3aWR0aCA9IDEuMikgKwogIGdlb21fZXJyb3JiYXIoYWVzKAogICAgeG1pbiA9IGVmZmVjdCAtIDEuOTYqc3RkLmVyciwKICAgIHhtYXggPSBlZmZlY3QgKyAxLjk2KnN0ZC5lcnIpLHdpZHRoID0gMC41LCBhbHBoYSA9IDAuNykgKwogIGdlb21fcG9pbnQoc2l6ZSA9IDEsY29sb3VyID0gImJsdWUiKSArIHRoZW1lX21pbmltYWwoKSArIAogIGxhYnMoeSA9ICJEaXN0cmljdCIsIHggPSAiUmFuZG9tIGVmZmVjdCIpICsgZ2d0aXRsZSgiUmFuZG9tIGVmZmVjdHMgZm9yIGVhY2ggcG9saWNlIGZvcmNlIikgKyAKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAibm9uZSIscGxvdC50aXRsZSA9IGVsZW1lbnRfdGV4dChzaXplID0gMTYsaGp1c3QgPSAwLjUpKSAKYGBgCgpgYGB7cn0KY2hlY2soYikKYGBgCgp3ZSBub3cgbmVlZCB0byBjb21wdXRlIHRoZSBiYW5kcyBpbiBhIGZ1bmN0aW9uYWwgY2FzZQoKYGBge3J9CnJlc190cmFpbiA8LSBjYmluZChkZl90cmFpbixwcmVkcyA9IHByZWRpY3QoZ2FtX3RyYWluLCBkZl90cmFpbiwgdHlwZSA9ICJyZXNwb25zZSIpKSAlPiUKICBhcnJhbmdlKHBvbGljZV9mb3JjZSxkYXRlKSAlPiUgbXV0YXRlKGFic2RpZmYgPSBhYnMobnVtYmVyX29mX2NyYXNoZXMtcHJlZHMpKQoKcmVzX2NhbCA8LSBjYmluZChkZl9jYWwscHJlZHMgPSBwcmVkaWN0KGdhbV90cmFpbiwgZGZfY2FsLCB0eXBlID0gInJlc3BvbnNlIikpICU+JQogIGFycmFuZ2UocG9saWNlX2ZvcmNlLGRhdGUpICU+JSBtdXRhdGUoYWJzZGlmZiA9IGFicyhudW1iZXJfb2ZfY3Jhc2hlcy1wcmVkcykpCgpyZXNfdGVzdCA8LSBjYmluZChkZl90ZXN0LHByZWRzID0gcHJlZGljdChnYW1fdHJhaW4sIGRmX3Rlc3QsIHR5cGUgPSAicmVzcG9uc2UiKSkgJT4lCiAgYXJyYW5nZShwb2xpY2VfZm9yY2UsZGF0ZSkgJT4lIG11dGF0ZShhYnNkaWZmID0gYWJzKG51bWJlcl9vZl9jcmFzaGVzLXByZWRzKSkKYGBgCgpgYGB7cn0Kbl9wb2wgPC0gbGVuZ3RoKGxldmVscyhjcmFzaGVzX3Blcl9kYXlfcG9saWNlJHBvbGljZV9mb3JjZSkpCgptID0gMTIKbCA9IDYKYWxwaGEgPSAwLjEKUyA9IG1hdHJpeChucm93ID0gbl9wb2wsbmNvbCA9IDM2NSkKayA9IHJlcCgwLG5fcG9sKQp4ID0gMTozNjUKbGluZV9pbnRlZ3JhbCA9IGZ1bmN0aW9uKHgsIHkpIHsKICBkeCA9IGRpZmYoeCkKICBlbmQgPSBsZW5ndGgoeSkKICBteSA9ICh5WzE6KGVuZCAtIDEpXSArIHlbMjplbmRdKSAvIDIKICBzdW0oZHggKm15KQp9IAoKY2VpbGluZygobSsxKSooMS1hbHBoYSkpCmBgYAoKd2UgdXNlIHRoZSBsYXJnZXN0IHZhbHVlIGZvciBnYW1tYQpIMSA9IEkxCgpjb21wdXRpbmcgdGhlIG1vZHVsYXRpb24gZnVuY3Rpb25zOgoKYGBge3J9CmZvcihpIGluIDE6bl9wb2wpewogIGFkX3RyYWluIDwtIHJlc190cmFpbiAlPiUgZmlsdGVyKHBvbGljZV9mb3JjZSA9PSBsZXZlbHMocG9saWNlX2ZvcmNlKVtpXSkgJT4lCiAgZHBseXI6OnNlbGVjdChhYnNkaWZmKQogIGFfbWF0IDwtIG1hdHJpeChhZF90cmFpbiRhYnNkaWZmLGJ5cm93ID0gRixucm93ID0gMzY1KQogIG5jbSA9IGFwcGx5KGFfbWF0LDEsbWF4KQogIGRlbiA9IGxpbmVfaW50ZWdyYWwoeCwgbmNtKQogIHMgPSBuY20vZGVuCiAgU1tpLF0gPSBzCn0KYGBgCgpjb21wdXRpbmcgdGhlIGs6CgpgYGB7cn0KY2VpbGluZygobCArIDEpKigxLWFscGhhKSkKYGBgCgp3ZSB1c2UgdGhlIGJpZ2dlc3QgdmFsdWUgZm9yIHRoZSBjaG9pY2Ugb2YgawoKYGBge3J9CmZvcihpIGluIDE6bl9wb2wpewogIGFkX2NhbCA8LSByZXNfY2FsICU+JSBmaWx0ZXIocG9saWNlX2ZvcmNlID09IGxldmVscyhwb2xpY2VfZm9yY2UpW2ldKSAlPiUKICBkcGx5cjo6c2VsZWN0KGFic2RpZmYpCiAgYV9tYXQgPC0gbWF0cml4KGFkX2NhbCRhYnNkaWZmLGJ5cm93ID0gRixucm93ID0gMzY1KQogIG5jcyA9IGFwcGx5KGFfbWF0L1NbaSxdLDEsbWF4KQogIGtbaV0gPSBtYXgobmNzKQp9CmBgYAoKYGBge3J9CmkgPC0gNTAKCmRmX2ZpbHRlcmVkIDwtIHJlc190ZXN0ICU+JSBmaWx0ZXIocG9saWNlX2ZvcmNlID09IGxldmVscyhwb2xpY2VfZm9yY2UpW2ldKSAlPiUKICBhcnJhbmdlKGRheV9vZl95ZWFyKQoKZGZfcGxvdCA8LSBkYXRhLmZyYW1lKG4gPSBkZl9maWx0ZXJlZCRudW1iZXJfb2ZfY3Jhc2hlcyxwcmVkID0gZGZfZmlsdGVyZWQkcHJlZHMsIAogICAgICAgICAgICAgICAgICAgICAgeW1heCA9IGRmX2ZpbHRlcmVkJHByZWRzICsga1tpXSpTW2ksXSwKICAgICAgICAgICAgICAgICAgICAgIHltaW4gPSBwbWF4KCBkZl9maWx0ZXJlZCRwcmVkcyAtIGtbaV0qU1tpLF0sIHJlcCgwLCAzNjUpKSx0ZCA9IDE6MzY1KQoKZGZfcGxvdCAlPiUgZ2dwbG90KCkgKyBnZW9tX3JpYmJvbihhZXMoeD10ZCx5ID0gcHJlZCx5bWluID0geW1pbiwgeW1heCA9IHltYXgpLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBmaWxsID0gImJsdWUiLCBhbHBoYSA9IDAuMjUpICsgCiAgZ2VvbV9saW5lKGFlcyh4PXRkLHkgPSBuKSxsaW5ld2lkdGg9MSxjb2xvdXIgPSAiYmxhY2siKSArIAogIGdlb21fbGluZShhZXMoeD10ZCx5ID0gcHJlZCksbGluZXdpZHRoPTEsY29sb3VyID0gImJsdWUiKSArIGxhYnMoeCA9ICJEYXkiLCB5ID0gIk51bWJlciBvZiBhY2NpZW50cyIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgdGl0bGUgPSBsZXZlbHMoZGZfdGVzdCRwb2xpY2VfZm9yY2UpW2ldKQoKYGBgCgp0aGlzIGlzIG5vdCBncmVhdDoKCgpgYGB7cn0KaSA8LSA2CgpkZl9maWx0ZXJlZCA8LSByZXNfdGVzdCAlPiUgZmlsdGVyKHBvbGljZV9mb3JjZSA9PSBsZXZlbHMocG9saWNlX2ZvcmNlKVtpXSkgJT4lCiAgYXJyYW5nZShkYXlfb2ZfeWVhcikKCmRmX3Bsb3QgPC0gZGF0YS5mcmFtZShuID0gZGZfZmlsdGVyZWQkbnVtYmVyX29mX2NyYXNoZXMscHJlZCA9IGRmX2ZpbHRlcmVkJHByZWRzLCAKICAgICAgICAgICAgICAgICAgICAgIHltYXggPSBkZl9maWx0ZXJlZCRwcmVkcyArIGtbaV0qU1tpLF0sCiAgICAgICAgICAgICAgICAgICAgICB5bWluID0gcG1heCggZGZfZmlsdGVyZWQkcHJlZHMgLSBrW2ldKlNbaSxdLCByZXAoMCwgMzY1KSksdGQgPSAxOjM2NSkKCmRmX3Bsb3QgJT4lIGdncGxvdCgpICsgZ2VvbV9yaWJib24oYWVzKHg9dGQseSA9IHByZWQseW1pbiA9IHltaW4sIHltYXggPSB5bWF4KSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZmlsbCA9ICJibHVlIiwgYWxwaGEgPSAwLjI1KSArIAogIGdlb21fbGluZShhZXMoeD10ZCx5ID0gbiksbGluZXdpZHRoPTEsY29sb3VyID0gImJsYWNrIikgKyAKICBnZW9tX2xpbmUoYWVzKHg9dGQseSA9IHByZWQpLGxpbmV3aWR0aD0xLGNvbG91ciA9ICJibHVlIikgKyBsYWJzKHggPSAiRGF5IiwgeSA9ICJOdW1iZXIgb2YgYWNjaWVudHMiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHRpdGxlID0gbGV2ZWxzKGRmX3Rlc3QkcG9saWNlX2ZvcmNlKVtpXSkKCmBgYAoKCmBgYHtyfQppIDwtIDMwCgpkZl9maWx0ZXJlZCA8LSByZXNfdGVzdCAlPiUgZmlsdGVyKHBvbGljZV9mb3JjZSA9PSBsZXZlbHMocG9saWNlX2ZvcmNlKVtpXSkgJT4lCiAgYXJyYW5nZShkYXlfb2ZfeWVhcikKCmRmX3Bsb3QgPC0gZGF0YS5mcmFtZShuID0gZGZfZmlsdGVyZWQkbnVtYmVyX29mX2NyYXNoZXMscHJlZCA9IGRmX2ZpbHRlcmVkJHByZWRzLCAKICAgICAgICAgICAgICAgICAgICAgIHltYXggPSBkZl9maWx0ZXJlZCRwcmVkcyArIGtbaV0qU1tpLF0sCiAgICAgICAgICAgICAgICAgICAgICB5bWluID0gcG1heCggZGZfZmlsdGVyZWQkcHJlZHMgLSBrW2ldKlNbaSxdLCByZXAoMCwgMzY1KSksdGQgPSAxOjM2NSkKCmRmX3Bsb3QgJT4lIGdncGxvdCgpICsgZ2VvbV9yaWJib24oYWVzKHg9dGQseSA9IHByZWQseW1pbiA9IHltaW4sIHltYXggPSB5bWF4KSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZmlsbCA9ICJibHVlIiwgYWxwaGEgPSAwLjI1KSArIAogIGdlb21fbGluZShhZXMoeD10ZCx5ID0gbiksbGluZXdpZHRoPTEsY29sb3VyID0gImJsYWNrIikgKyAKICBnZW9tX2xpbmUoYWVzKHg9dGQseSA9IHByZWQpLGxpbmV3aWR0aD0xLGNvbG91ciA9ICJibHVlIikgKyBsYWJzKHggPSAiRGF5IiwgeSA9ICJOdW1iZXIgb2YgYWNjaWVudHMiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHRpdGxlID0gbGV2ZWxzKGRmX3Rlc3QkcG9saWNlX2ZvcmNlKVtpXSkKCmBgYAoKdGhpcyBhcmUgdG9vIGJpZywgaXQgaXMgZHVlIHRvIHRoZSBmYWN0IHRoYXQgd2UgYXJlIHVzaW5nIGEgbWl4ZCBlZmZlY3RzIG1vZGVsLCAKc28gbWV0cm9wb2xpdGFuIHBvbGljZSB3aWxsIGNvbnRhbWluYXRlIGFsbCBvZiB0aGUgZGF0YS4KCndlIGNhbiBmaXggdGhpcyBieSBtYWtpbmcgYSBnYW0gZm9yIGVhY2ggcG9saWNlIGZvcmNlOgoKYGBge3J9CmdhbXMgPC0gbGlzdCgpCgpTID0gbWF0cml4KG5yb3cgPSBuX3BvbCxuY29sID0gMzY1KQprID0gcmVwKDAsbl9wb2wpCgpmb3IoaSBpbiAxOm5fcG9sKXsKICAKICBkZl90cmFpbl9taWQgPC0gZGZfdHJhaW4gJT4lIGZpbHRlcihwb2xpY2VfZm9yY2UgPT0gbGV2ZWxzKHBvbGljZV9mb3JjZSlbaV0pCiAgZGZfY2FsX21pZCA8LSBkZl9jYWwgJT4lIGZpbHRlcihwb2xpY2VfZm9yY2UgPT0gbGV2ZWxzKHBvbGljZV9mb3JjZSlbaV0pCiAgCiAgZ2FtIDwtIGJhbShudW1iZXJfb2ZfY3Jhc2hlcyB+IGRheSArIHllYXIgKwogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcyhkYXlfb2ZfeWVhcixrID0gNTMsIGJzID0gImNyIiksIAogICAgICAgICAgICAgICAgZGF0YSA9IGRmX3RyYWluX21pZCwgZmFtaWx5PXF1YXNpcG9pc3NvbigpLCBtZXRob2Q9J1JFTUwnKQogIAogIHJlc190cmFpbl9taWQgPC0gY2JpbmQoZGZfdHJhaW5fbWlkLHByZWRzID0gcHJlZGljdChnYW0sIGRmX3RyYWluX21pZCwgdHlwZSA9ICJyZXNwb25zZSIpKSAlPiUKICBhcnJhbmdlKGRhdGUpICU+JSBtdXRhdGUoYWJzZGlmZiA9IGFicyhudW1iZXJfb2ZfY3Jhc2hlcy1wcmVkcykpCgogIHJlc19jYWxfbWlkIDwtIGNiaW5kKGRmX2NhbF9taWQscHJlZHMgPSBwcmVkaWN0KGdhbSwgZGZfY2FsX21pZCwgdHlwZSA9ICJyZXNwb25zZSIpKSAlPiUKICBhcnJhbmdlKGRhdGUpICU+JSBtdXRhdGUoYWJzZGlmZiA9IGFicyhudW1iZXJfb2ZfY3Jhc2hlcy1wcmVkcykpCiAgCiAgCiAgYWRfdHJhaW4gPC0gcmVzX3RyYWluX21pZCAlPiUgZmlsdGVyKHBvbGljZV9mb3JjZSA9PSBsZXZlbHMocG9saWNlX2ZvcmNlKVtpXSkgJT4lCiAgZHBseXI6OnNlbGVjdChhYnNkaWZmKQogIGFfbWF0IDwtIG1hdHJpeChhZF90cmFpbiRhYnNkaWZmLGJ5cm93ID0gRixucm93ID0gMzY1KQogIG5jbSA9IGFwcGx5KGFfbWF0LDEsbWF4KQogIGRlbiA9IGxpbmVfaW50ZWdyYWwoeCwgbmNtKQogIHMgPSBuY20vZGVuCiAgU1tpLF0gPSBzCiAgYWRfY2FsIDwtIHJlc19jYWxfbWlkICU+JSBmaWx0ZXIocG9saWNlX2ZvcmNlID09IGxldmVscyhwb2xpY2VfZm9yY2UpW2ldKSAlPiUKICBkcGx5cjo6c2VsZWN0KGFic2RpZmYpCiAgYV9tYXQgPC0gbWF0cml4KGFkX2NhbCRhYnNkaWZmLGJ5cm93ID0gRixucm93ID0gMzY1KQogIG5jcyA9IGFwcGx5KGFfbWF0L1NbaSxdLDEsbWF4KQogIGtbaV0gPSBtYXgobmNzKQogIAogIGdhbXNbW2ldXSA8LSBnYW0KfQpgYGAKCm1pZGxhbmRzCgpgYGB7cn0KaSA8LSA1MApkZl90ZXN0X2ZpbHRlcmVkIDwtIGRmX3Rlc3QgJT4lIGZpbHRlcihwb2xpY2VfZm9yY2UgPT0gbGV2ZWxzKHBvbGljZV9mb3JjZSlbaV0pCnByZWRzIDwtIHByZWRpY3QoZ2Ftc1tbaV1dLCBkZl90ZXN0X2ZpbHRlcmVkLCB0eXBlID0gInJlc3BvbnNlIikKZGZfdGVzdF9maWx0ZXJlZCAgPC0gY2JpbmQoZGZfdGVzdF9maWx0ZXJlZCxwcmVkcykgJT4lIGFycmFuZ2UoZGF5X29mX3llYXIpCmBgYAoKYGBge3J9CmRmX3Bsb3QgPC0gZGF0YS5mcmFtZShuID0gZGZfdGVzdF9maWx0ZXJlZCRudW1iZXJfb2ZfY3Jhc2hlcyxwcmVkID0gZGZfdGVzdF9maWx0ZXJlZCRwcmVkcywgCiAgICAgICAgICAgICAgICAgICAgICB5bWF4ID0gZGZfdGVzdF9maWx0ZXJlZCRwcmVkcyArIGtbaV0qU1tpLF0sCiAgICAgICAgICAgICAgICAgICAgICB5bWluID0gcG1heCggZGZfdGVzdF9maWx0ZXJlZCRwcmVkcyAtIGtbaV0qU1tpLF0sIHJlcCgwLCAzNjUpKSx0ZCA9IDE6MzY1KQoKZGZfcGxvdCAlPiUgZ2dwbG90KCkgKyBnZW9tX3JpYmJvbihhZXMoeD10ZCx5ID0gcHJlZCx5bWluID0geW1pbiwgeW1heCA9IHltYXgpLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBmaWxsID0gImJsdWUiLCBhbHBoYSA9IDAuMjUpICsgCiAgZ2VvbV9saW5lKGFlcyh4PXRkLHkgPSBuKSxsaW5ld2lkdGg9MSxjb2xvdXIgPSAiYmxhY2siKSArIAogIGdlb21fbGluZShhZXMoeD10ZCx5ID0gcHJlZCksbGluZXdpZHRoPTEsY29sb3VyID0gImJsdWUiKSArIGxhYnMoeCA9ICJEYXkiLCB5ID0gIk51bWJlciBvZiBhY2NpZW50cyIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgdGl0bGUgPSBsZXZlbHMoZGZfdGVzdCRwb2xpY2VfZm9yY2UpW2ldKQpgYGAKCkxvbmRvbjoKCmBgYHtyfQppIDwtIDMwCmRmX3Rlc3RfZmlsdGVyZWQgPC0gZGZfdGVzdCAlPiUgZmlsdGVyKHBvbGljZV9mb3JjZSA9PSBsZXZlbHMocG9saWNlX2ZvcmNlKVtpXSkKcHJlZHMgPC0gcHJlZGljdChnYW1zW1tpXV0sIGRmX3Rlc3RfZmlsdGVyZWQsIHR5cGUgPSAicmVzcG9uc2UiKQpkZl90ZXN0X2ZpbHRlcmVkICA8LSBjYmluZChkZl90ZXN0X2ZpbHRlcmVkLHByZWRzKSAlPiUgYXJyYW5nZShkYXlfb2ZfeWVhcikKYGBgCgpgYGB7cn0KZGZfcGxvdCA8LSBkYXRhLmZyYW1lKG4gPSBkZl90ZXN0X2ZpbHRlcmVkJG51bWJlcl9vZl9jcmFzaGVzLHByZWQgPSBkZl90ZXN0X2ZpbHRlcmVkJHByZWRzLCAKICAgICAgICAgICAgICAgICAgICAgIHltYXggPSBkZl90ZXN0X2ZpbHRlcmVkJHByZWRzICsga1tpXSpTW2ksXSwKICAgICAgICAgICAgICAgICAgICAgIHltaW4gPSBwbWF4KCBkZl90ZXN0X2ZpbHRlcmVkJHByZWRzIC0ga1tpXSpTW2ksXSwgcmVwKDAsIDM2NSkpLHRkID0gMTozNjUpCgpkZl9wbG90ICU+JSBnZ3Bsb3QoKSArIGdlb21fcmliYm9uKGFlcyh4PXRkLHkgPSBwcmVkLHltaW4gPSB5bWluLCB5bWF4ID0geW1heCksIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZpbGwgPSAiYmx1ZSIsIGFscGhhID0gMC4yNSkgKyAKICBnZW9tX2xpbmUoYWVzKHg9dGQseSA9IG4pLGxpbmV3aWR0aD0xLGNvbG91ciA9ICJibGFjayIpICsgCiAgZ2VvbV9saW5lKGFlcyh4PXRkLHkgPSBwcmVkKSxsaW5ld2lkdGg9MSxjb2xvdXIgPSAiYmx1ZSIpICsgbGFicyh4ID0gIkRheSIsIHkgPSAiTnVtYmVyIG9mIGFjY2llbnRzIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB0aXRsZSA9IGxldmVscyhkZl90ZXN0JHBvbGljZV9mb3JjZSlbaV0pCmBgYAoKdGhlIG1vZGVsIGlzIGJldHRlciBidXQgc3RpbGwgbm8gZ3JlYXQKCndlIHByb2JhYmx5IG5lZWQgbW9yZSBkYXRhLgoKdXNpbmcgYWxsIHRoZSBhdmFpbGFibGUgZGF0YToKCmBgYHtyfQpmdWxsX2NvbGxpc2lvaW5fZGF0YSA8LSByZWFkcjo6cmVhZF9jc3YoZmlsZT0iZGZ0LXJvYWQtY2FzdWFsdHktc3RhdGlzdGljcy1jb2xsaXNpb24tMTk3OS1sYXRlc3QtcHVibGlzaGVkLXllYXIuY3N2IikgCgpzcGVjKGZ1bGxfY29sbGlzaW9pbl9kYXRhKQpnbGltcHNlKGZ1bGxfY29sbGlzaW9pbl9kYXRhKQoKZnVsbF9jb2xsaXNpb25zX2FsbCA8LSBzdGF0czE5Ojpmb3JtYXRfY29sbGlzaW9ucyhmdWxsX2NvbGxpc2lvaW5fZGF0YSkKCmdsaW1wc2UoZnVsbF9jb2xsaXNpb25zX2FsbCkKYGBgCgpgYGB7cn0KYWxsX2RhdGVzIDwtIHVuaXF1ZShmdWxsX2NvbGxpc2lvbnNfYWxsJGRhdGUpCmFsbF9wb2xpY2VfZm9yY2VzIDwtIHVuaXF1ZShmdWxsX2NvbGxpc2lvbnNfYWxsJHBvbGljZV9mb3JjZSkKCmZ1bGxfY29tYmluYXRpb25zIDwtIGV4cGFuZC5ncmlkKGRhdGUgPSBhbGxfZGF0ZXMsIHBvbGljZV9mb3JjZSA9IGFsbF9wb2xpY2VfZm9yY2VzKQoKY3Jhc2hlc19wZXJfZGF5X3BvbGljZSA8LSBmdWxsX2NvbWJpbmF0aW9ucyAlPiUKICBsZWZ0X2pvaW4oZnVsbF9jb2xsaXNpb25zX2FsbCAlPiUgZ3JvdXBfYnkoZGF0ZSwgcG9saWNlX2ZvcmNlKSAlPiUgCiAgICAgICAgICAgICAgc3VtbWFyaXplKG51bWJlcl9vZl9jcmFzaGVzID0gbigpKSwKICAgICAgICAgICAgYnkgPSBjKCJkYXRlIiwgInBvbGljZV9mb3JjZSIpKSAlPiUKICByZXBsYWNlX25hKGxpc3QobnVtYmVyX29mX2NyYXNoZXMgPSAwKSkKCmdsaW1wc2UoY3Jhc2hlc19wZXJfZGF5X3BvbGljZSkKYGBgCgpgYGB7cn0KY3Jhc2hlc19wZXJfZGF5X3BvbGljZSA8LSBjcmFzaGVzX3Blcl9kYXlfcG9saWNlICU+JSBtdXRhdGUoZGF5X29mX3llYXIgPSB5ZGF5KGRhdGUpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBkYXkgPSBmYWN0b3Iod2RheShkYXRlKSkseWVhciA9IHllYXIoZGF0ZSkpCmBgYAoKcmVtb3ZpbmcgdGhlIDI5IG9mIGZlYnJ1YXJ5IGZvciBkaW1lbnNpb2FsaXR5IHByb2JsZW1zIGluIHRoZSB2ZWN0b3JzIGJlbG93CgpgYGB7cn0KY3Jhc2hlc19wZXJfZGF5X3BvbGljZSA8LSBjcmFzaGVzX3Blcl9kYXlfcG9saWNlICU+JSBmaWx0ZXIoIShtb250aChkYXRlKT09MiAmIGRheShkYXRlKT09MjkpKQpgYGAKCmBgYHtyfQp0cmFpbl95ZWFycyA8LSBjZWlsaW5nKHNlcSgxOTc5LDIwMjEsYnkgPSAxLjUpKQpjYWxfeWVhcnMgPC0gc2VxKDE5NzksMjAyMSwzKQoKZGZfdHJhaW4gPC0gY3Jhc2hlc19wZXJfZGF5X3BvbGljZSAlPiUgZmlsdGVyKHllYXIgJWluJSB0cmFpbl95ZWFycykgCmRmX2NhbCA8LSBjcmFzaGVzX3Blcl9kYXlfcG9saWNlICU+JSBmaWx0ZXIoeWVhciAlaW4lIGNhbF95ZWFycykgCmRmX3Rlc3QgPC0gY3Jhc2hlc19wZXJfZGF5X3BvbGljZSAlPiUgZmlsdGVyKHllYXIgPT0gMjAyMikKYGBgCgpzdGFydGluZyB3aXRoIGEgZ2FtIGFuZCBtaXhlZCBlZmZlY3RzCgpgYGB7cn0KZ2FtX3RyYWluIDwtIGJhbShudW1iZXJfb2ZfY3Jhc2hlcyB+IGRheSArIHllYXIgKwogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcyhkYXlfb2ZfeWVhcixrID0gNTMsIGJzID0gImNyIikgKyAKICAgICAgICAgICAgICAgICAgICAgICBzKHBvbGljZV9mb3JjZSwgYnMgPSAicmUiKSwgCiAgICAgICAgICAgICAgICBkYXRhID0gZGZfdHJhaW4sIGZhbWlseT1xdWFzaXBvaXNzb24oKSwgbWV0aG9kPSdSRU1MJykKYGBgCgpgYGB7cn0Kc3VtbWFyeShnYW1fdHJhaW4pCmBgYAoKYGBge3J9CmIgPC0gZ2V0Vml6KGdhbV90cmFpbikKYGBgCgpgYGB7cn0KcHJpbnQocGxvdChiLCBhbGxUZXJtcyA9IFQpLCBwYWdlcyA9IDEpCmBgYAoKMSBpcyBzdW5kYXkgaGVyZQoKYGBge3J9Cm5ld2QgPC0gd2l0aChkZl90cmFpbiwgZGF0YS5mcmFtZShkYXlfb2ZfeWVhciA9IDEyMCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGRheSA9IDEsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5ZWFyID0gMjAyMiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHBvbGljZV9mb3JjZSA9IGxldmVscyhwb2xpY2VfZm9yY2UpKSkKCnAgPC0gcHJlZGljdChnYW1fdHJhaW4sIG5ld2QsIHR5cGUgPSAidGVybXMiLCBzZS5maXQgPSBUUlVFKQpyZSA8LSBwW1siZml0Il1dWyAsICJzKHBvbGljZV9mb3JjZSkiXQpzZSA8LSBwW1sic2UuZml0Il1dWyAsICJzKHBvbGljZV9mb3JjZSkiXQoKZGF0YSA8LSBkYXRhLmZyYW1lKHBvbGljZV9mb3JjZSA9IGxldmVscyhkZl90cmFpbiRwb2xpY2VfZm9yY2UpLCBlZmZlY3QgPSByZSxzdGQuZXJyID0gc2UpCmBgYAoKYGBge3J9CmRhdGEgJT4lIApnZ3Bsb3QoYWVzKGVmZmVjdCwgZmN0X3Jlb3JkZXIocG9saWNlX2ZvcmNlLCBlZmZlY3QpKSkgKwogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IDAsIGNvbG9yID0gImdyYXk1MCIsIGx0eSA9IDIsIGxpbmV3aWR0aCA9IDEuMikgKwogIGdlb21fZXJyb3JiYXIoYWVzKAogICAgeG1pbiA9IGVmZmVjdCAtIDEuOTYqc3RkLmVyciwKICAgIHhtYXggPSBlZmZlY3QgKyAxLjk2KnN0ZC5lcnIpLHdpZHRoID0gMC41LCBhbHBoYSA9IDAuNykgKwogIGdlb21fcG9pbnQoc2l6ZSA9IDEsY29sb3VyID0gImJsdWUiKSArIHRoZW1lX21pbmltYWwoKSArIAogIGxhYnMoeSA9ICJEaXN0cmljdCIsIHggPSAiUmFuZG9tIGVmZmVjdCIpICsgZ2d0aXRsZSgiUmFuZG9tIGVmZmVjdHMgZm9yIGVhY2ggcG9saWNlIGZvcmNlIikgKyAKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAibm9uZSIscGxvdC50aXRsZSA9IGVsZW1lbnRfdGV4dChzaXplID0gMTYsaGp1c3QgPSAwLjUpKSAKYGBgCgpgYGB7cn0KY2hlY2soYikKYGBgCgp3ZSBub3cgbmVlZCB0byBjb21wdXRlIHRoZSBiYW5kcyBpbiBhIGZ1bmN0aW9uYWwgY2FzZQoKYGBge3J9CnJlc190cmFpbiA8LSBjYmluZChkZl90cmFpbixwcmVkcyA9IHByZWRpY3QoZ2FtX3RyYWluLCBkZl90cmFpbiwgdHlwZSA9ICJyZXNwb25zZSIpKSAlPiUKICBhcnJhbmdlKHBvbGljZV9mb3JjZSxkYXRlKSAlPiUgbXV0YXRlKGFic2RpZmYgPSBhYnMobnVtYmVyX29mX2NyYXNoZXMtcHJlZHMpKQoKcmVzX2NhbCA8LSBjYmluZChkZl9jYWwscHJlZHMgPSBwcmVkaWN0KGdhbV90cmFpbiwgZGZfY2FsLCB0eXBlID0gInJlc3BvbnNlIikpICU+JQogIGFycmFuZ2UocG9saWNlX2ZvcmNlLGRhdGUpICU+JSBtdXRhdGUoYWJzZGlmZiA9IGFicyhudW1iZXJfb2ZfY3Jhc2hlcy1wcmVkcykpCgpyZXNfdGVzdCA8LSBjYmluZChkZl90ZXN0LHByZWRzID0gcHJlZGljdChnYW1fdHJhaW4sIGRmX3Rlc3QsIHR5cGUgPSAicmVzcG9uc2UiKSkgJT4lCiAgYXJyYW5nZShwb2xpY2VfZm9yY2UsZGF0ZSkgJT4lIG11dGF0ZShhYnNkaWZmID0gYWJzKG51bWJlcl9vZl9jcmFzaGVzLXByZWRzKSkKYGBgCgpgYGB7cn0Kbl9wb2wgPC0gbGVuZ3RoKGxldmVscyhjcmFzaGVzX3Blcl9kYXlfcG9saWNlJHBvbGljZV9mb3JjZSkpCgptID0gbGVuZ3RoKHRyYWluX3llYXJzKQpsID0gbGVuZ3RoKGNhbF95ZWFycykKYWxwaGEgPSAwLjEKUyA9IG1hdHJpeChucm93ID0gbl9wb2wsbmNvbCA9IDM2NSkKayA9IHJlcCgwLG5fcG9sKQp4ID0gMTozNjUKbGluZV9pbnRlZ3JhbCA9IGZ1bmN0aW9uKHgsIHkpIHsKICBkeCA9IGRpZmYoeCkKICBlbmQgPSBsZW5ndGgoeSkKICBteSA9ICh5WzE6KGVuZCAtIDEpXSArIHlbMjplbmRdKSAvIDIKICBzdW0oZHggKm15KQp9IAoKY2VpbGluZygobSsxKSooMS1hbHBoYSkpCmBgYAoKSDEgIT0gSTEgbm93Cgpjb21wdXRpbmcgdGhlIG1vZHVsYXRpb24gZnVuY3Rpb25zOgoKYGBge3J9CmZvcihpIGluIDE6bl9wb2wpewogIGFkX3RyYWluIDwtIHJlc190cmFpbiAlPiUgZmlsdGVyKHBvbGljZV9mb3JjZSA9PSBsZXZlbHMocG9saWNlX2ZvcmNlKVtpXSkgJT4lCiAgZHBseXI6OnNlbGVjdChhYnNkaWZmKQogIGFfbWF0IDwtIG1hdHJpeChhZF90cmFpbiRhYnNkaWZmLGJ5cm93ID0gRixucm93ID0gMzY1KQogIHNldGNmciA8LSBhcHBseShhX21hdCwyLG1heCkKICBnYW1tYSA8LSBzb3J0KHNldGNmcilbY2VpbGluZygobSsxKSooMS1hbHBoYSkpXQogIGFfbWF0X2ZpbHQgPC0gYV9tYXRbLHNldGNmcjw9Z2FtbWFdCiAgbmNtID0gYXBwbHkoYV9tYXRfZmlsdCwxLG1heCkKICBuY21fc29ydCA8LSBzb3J0KG5jbSkKICBkZW4gPSBsaW5lX2ludGVncmFsKHgsIG5jbSkKICBzID0gbmNtL2RlbgogIFNbaSxdID0gcwp9CmBgYAoKY29tcHV0aW5nIHRoZSBrOgoKYGBge3J9CmNlaWxpbmcoKGwgKyAxKSooMS1hbHBoYSkpCmBgYAoKd2UgdXNlIHRoZSBiaWdnZXN0IHZhbHVlIGZvciB0aGUgY2hvaWNlIG9mIGsKCmBgYHtyfQpmb3IoaSBpbiAxOm5fcG9sKXsKICBhZF9jYWwgPC0gcmVzX2NhbCAlPiUgZmlsdGVyKHBvbGljZV9mb3JjZSA9PSBsZXZlbHMocG9saWNlX2ZvcmNlKVtpXSkgJT4lCiAgZHBseXI6OnNlbGVjdChhYnNkaWZmKQogIGFfbWF0IDwtIG1hdHJpeChhZF9jYWwkYWJzZGlmZixieXJvdyA9IEYsbnJvdyA9IDM2NSkKICBuY3MgPSBhcHBseShhX21hdC9TW2ksXSwxLG1heCkKICBrW2ldID0gbWF4KG5jcykKfQpgYGAKCmBgYHtyfQppIDwtIDUwCgpkZl9maWx0ZXJlZCA8LSByZXNfdGVzdCAlPiUgZmlsdGVyKHBvbGljZV9mb3JjZSA9PSBsZXZlbHMocG9saWNlX2ZvcmNlKVtpXSkgJT4lCiAgYXJyYW5nZShkYXlfb2ZfeWVhcikKCmRmX3Bsb3QgPC0gZGF0YS5mcmFtZShuID0gZGZfZmlsdGVyZWQkbnVtYmVyX29mX2NyYXNoZXMscHJlZCA9IGRmX2ZpbHRlcmVkJHByZWRzLCAKICAgICAgICAgICAgICAgICAgICAgIHltYXggPSBkZl9maWx0ZXJlZCRwcmVkcyArIGtbaV0qU1tpLF0sCiAgICAgICAgICAgICAgICAgICAgICB5bWluID0gcG1heCggZGZfZmlsdGVyZWQkcHJlZHMgLSBrW2ldKlNbaSxdLCByZXAoMCwgMzY1KSksdGQgPSAxOjM2NSkKCmRmX3Bsb3QgJT4lIGdncGxvdCgpICsgZ2VvbV9yaWJib24oYWVzKHg9dGQseSA9IHByZWQseW1pbiA9IHltaW4sIHltYXggPSB5bWF4KSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZmlsbCA9ICJibHVlIiwgYWxwaGEgPSAwLjI1KSArIAogIGdlb21fbGluZShhZXMoeD10ZCx5ID0gbiksbGluZXdpZHRoPTEsY29sb3VyID0gImJsYWNrIikgKyAKICBnZW9tX2xpbmUoYWVzKHg9dGQseSA9IHByZWQpLGxpbmV3aWR0aD0xLGNvbG91ciA9ICJibHVlIikgKyBsYWJzKHggPSAiRGF5IiwgeSA9ICJOdW1iZXIgb2YgYWNjaWVudHMiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHRpdGxlID0gbGV2ZWxzKGRmX3Rlc3QkcG9saWNlX2ZvcmNlKVtpXSkKCmBgYAoKdGhpcyBpcyBub3QgZ3JlYXQ6CgoKYGBge3J9CmkgPC0gNgoKZGZfZmlsdGVyZWQgPC0gcmVzX3Rlc3QgJT4lIGZpbHRlcihwb2xpY2VfZm9yY2UgPT0gbGV2ZWxzKHBvbGljZV9mb3JjZSlbaV0pICU+JQogIGFycmFuZ2UoZGF5X29mX3llYXIpCgpkZl9wbG90IDwtIGRhdGEuZnJhbWUobiA9IGRmX2ZpbHRlcmVkJG51bWJlcl9vZl9jcmFzaGVzLHByZWQgPSBkZl9maWx0ZXJlZCRwcmVkcywgCiAgICAgICAgICAgICAgICAgICAgICB5bWF4ID0gZGZfZmlsdGVyZWQkcHJlZHMgKyBrW2ldKlNbaSxdLAogICAgICAgICAgICAgICAgICAgICAgeW1pbiA9IHBtYXgoIGRmX2ZpbHRlcmVkJHByZWRzIC0ga1tpXSpTW2ksXSwgcmVwKDAsIDM2NSkpLHRkID0gMTozNjUpCgpkZl9wbG90ICU+JSBnZ3Bsb3QoKSArIGdlb21fcmliYm9uKGFlcyh4PXRkLHkgPSBwcmVkLHltaW4gPSB5bWluLCB5bWF4ID0geW1heCksIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZpbGwgPSAiYmx1ZSIsIGFscGhhID0gMC4yNSkgKyAKICBnZW9tX2xpbmUoYWVzKHg9dGQseSA9IG4pLGxpbmV3aWR0aD0xLGNvbG91ciA9ICJibGFjayIpICsgCiAgZ2VvbV9saW5lKGFlcyh4PXRkLHkgPSBwcmVkKSxsaW5ld2lkdGg9MSxjb2xvdXIgPSAiYmx1ZSIpICsgbGFicyh4ID0gIkRheSIsIHkgPSAiTnVtYmVyIG9mIGFjY2llbnRzIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB0aXRsZSA9IGxldmVscyhkZl90ZXN0JHBvbGljZV9mb3JjZSlbaV0pCgpgYGAKCgpgYGB7cn0KaSA8LSAzMAoKZGZfZmlsdGVyZWQgPC0gcmVzX3Rlc3QgJT4lIGZpbHRlcihwb2xpY2VfZm9yY2UgPT0gbGV2ZWxzKHBvbGljZV9mb3JjZSlbaV0pICU+JQogIGFycmFuZ2UoZGF5X29mX3llYXIpCgpkZl9wbG90IDwtIGRhdGEuZnJhbWUobiA9IGRmX2ZpbHRlcmVkJG51bWJlcl9vZl9jcmFzaGVzLHByZWQgPSBkZl9maWx0ZXJlZCRwcmVkcywgCiAgICAgICAgICAgICAgICAgICAgICB5bWF4ID0gZGZfZmlsdGVyZWQkcHJlZHMgKyBrW2ldKlNbaSxdLAogICAgICAgICAgICAgICAgICAgICAgeW1pbiA9IHBtYXgoIGRmX2ZpbHRlcmVkJHByZWRzIC0ga1tpXSpTW2ksXSwgcmVwKDAsIDM2NSkpLHRkID0gMTozNjUpCgpkZl9wbG90ICU+JSBnZ3Bsb3QoKSArIGdlb21fcmliYm9uKGFlcyh4PXRkLHkgPSBwcmVkLHltaW4gPSB5bWluLCB5bWF4ID0geW1heCksIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZpbGwgPSAiYmx1ZSIsIGFscGhhID0gMC4yNSkgKyAKICBnZW9tX2xpbmUoYWVzKHg9dGQseSA9IG4pLGxpbmV3aWR0aD0xLGNvbG91ciA9ICJibGFjayIpICsgCiAgZ2VvbV9saW5lKGFlcyh4PXRkLHkgPSBwcmVkKSxsaW5ld2lkdGg9MSxjb2xvdXIgPSAiYmx1ZSIpICsgbGFicyh4ID0gIkRheSIsIHkgPSAiTnVtYmVyIG9mIGFjY2llbnRzIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB0aXRsZSA9IGxldmVscyhkZl90ZXN0JHBvbGljZV9mb3JjZSlbaV0pCgpgYGAKCnRoaXMgaXMgZXZlIHdvcnNlIHRoYSBiZWZvcmU6Cgpkb2luZyBpdCBmb3IgYSBzaW5nbGUgcG9saWNlIGZvcmNlIGF0IGEgdGltZToKCmBgYHtyfQpnYW1zIDwtIGxpc3QoKQoKUyA9IG1hdHJpeChucm93ID0gbl9wb2wsbmNvbCA9IDM2NSkKayA9IHJlcCgwLG5fcG9sKQoKZm9yKGkgaW4gMTpuX3BvbCl7CiAgCiAgZGZfdHJhaW5fbWlkIDwtIGRmX3RyYWluICU+JSBmaWx0ZXIocG9saWNlX2ZvcmNlID09IGxldmVscyhwb2xpY2VfZm9yY2UpW2ldKQogIGRmX2NhbF9taWQgPC0gZGZfY2FsICU+JSBmaWx0ZXIocG9saWNlX2ZvcmNlID09IGxldmVscyhwb2xpY2VfZm9yY2UpW2ldKQogIAogIGdhbSA8LSBiYW0obnVtYmVyX29mX2NyYXNoZXMgfiBkYXkgKyB5ZWFyICsKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHMoZGF5X29mX3llYXIsayA9IDUzLCBicyA9ICJjciIpLCAKICAgICAgICAgICAgICAgIGRhdGEgPSBkZl90cmFpbl9taWQsIGZhbWlseT1xdWFzaXBvaXNzb24oKSwgbWV0aG9kPSdSRU1MJykKICAKICByZXNfdHJhaW5fbWlkIDwtIGNiaW5kKGRmX3RyYWluX21pZCxwcmVkcyA9IHByZWRpY3QoZ2FtLCBkZl90cmFpbl9taWQsIHR5cGUgPSAicmVzcG9uc2UiKSkgJT4lCiAgYXJyYW5nZShkYXRlKSAlPiUgbXV0YXRlKGFic2RpZmYgPSBhYnMobnVtYmVyX29mX2NyYXNoZXMtcHJlZHMpKQoKICByZXNfY2FsX21pZCA8LSBjYmluZChkZl9jYWxfbWlkLHByZWRzID0gcHJlZGljdChnYW0sIGRmX2NhbF9taWQsIHR5cGUgPSAicmVzcG9uc2UiKSkgJT4lCiAgYXJyYW5nZShkYXRlKSAlPiUgbXV0YXRlKGFic2RpZmYgPSBhYnMobnVtYmVyX29mX2NyYXNoZXMtcHJlZHMpKQogIAogIGFkX3RyYWluIDwtIHJlc190cmFpbl9taWQgJT4lIGZpbHRlcihwb2xpY2VfZm9yY2UgPT0gbGV2ZWxzKHBvbGljZV9mb3JjZSlbaV0pICU+JQogIGRwbHlyOjpzZWxlY3QoYWJzZGlmZikKICBhX21hdCA8LSBtYXRyaXgoYWRfdHJhaW4kYWJzZGlmZixieXJvdyA9IEYsbnJvdyA9IDM2NSkKICBzZXRjZnIgPC0gYXBwbHkoYV9tYXQsMixtYXgpCiAgZ2FtbWEgPC0gc29ydChzZXRjZnIpW2NlaWxpbmcoKG0rMSkqKDEtYWxwaGEpKV0KICBhX21hdF9maWx0IDwtIGFfbWF0WyxzZXRjZnI8PWdhbW1hXQogIG5jbSA9IGFwcGx5KGFfbWF0X2ZpbHQsMSxtYXgpCiAgbmNtX3NvcnQgPC0gc29ydChuY20pCiAgZGVuID0gbGluZV9pbnRlZ3JhbCh4LCBuY20pCiAgcyA9IG5jbS9kZW4KICBTW2ksXSA9IHMKICAKICBhZF9jYWwgPC0gcmVzX2NhbF9taWQgJT4lIGZpbHRlcihwb2xpY2VfZm9yY2UgPT0gbGV2ZWxzKHBvbGljZV9mb3JjZSlbaV0pICU+JQogIGRwbHlyOjpzZWxlY3QoYWJzZGlmZikKICBhX21hdCA8LSBtYXRyaXgoYWRfY2FsJGFic2RpZmYsYnlyb3cgPSBGLG5yb3cgPSAzNjUpCiAgbmNzID0gYXBwbHkoYV9tYXQvU1tpLF0sMSxtYXgpCiAga1tpXSA9IG1heChuY3MpCiAgCiAgZ2Ftc1tbaV1dIDwtIGdhbQp9CmBgYAoKbWlkbGFuZHMKCmBgYHtyfQppIDwtIDUwCmRmX3Rlc3RfZmlsdGVyZWQgPC0gZGZfdGVzdCAlPiUgZmlsdGVyKHBvbGljZV9mb3JjZSA9PSBsZXZlbHMocG9saWNlX2ZvcmNlKVtpXSkKcHJlZHMgPC0gcHJlZGljdChnYW1zW1tpXV0sIGRmX3Rlc3RfZmlsdGVyZWQsIHR5cGUgPSAicmVzcG9uc2UiKQpkZl90ZXN0X2ZpbHRlcmVkICA8LSBjYmluZChkZl90ZXN0X2ZpbHRlcmVkLHByZWRzKSAlPiUgYXJyYW5nZShkYXlfb2ZfeWVhcikKYGBgCgpgYGB7cn0KZGZfcGxvdCA8LSBkYXRhLmZyYW1lKG4gPSBkZl90ZXN0X2ZpbHRlcmVkJG51bWJlcl9vZl9jcmFzaGVzLHByZWQgPSBkZl90ZXN0X2ZpbHRlcmVkJHByZWRzLCAKICAgICAgICAgICAgICAgICAgICAgIHltYXggPSBkZl90ZXN0X2ZpbHRlcmVkJHByZWRzICsga1tpXSpTW2ksXSwKICAgICAgICAgICAgICAgICAgICAgIHltaW4gPSBwbWF4KCBkZl90ZXN0X2ZpbHRlcmVkJHByZWRzIC0ga1tpXSpTW2ksXSwgcmVwKDAsIDM2NSkpLHRkID0gMTozNjUpCgpkZl9wbG90ICU+JSBnZ3Bsb3QoKSArIGdlb21fcmliYm9uKGFlcyh4PXRkLHkgPSBwcmVkLHltaW4gPSB5bWluLCB5bWF4ID0geW1heCksIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZpbGwgPSAiYmx1ZSIsIGFscGhhID0gMC4yNSkgKyAKICBnZW9tX2xpbmUoYWVzKHg9dGQseSA9IG4pLGxpbmV3aWR0aD0xLGNvbG91ciA9ICJibGFjayIpICsgCiAgZ2VvbV9saW5lKGFlcyh4PXRkLHkgPSBwcmVkKSxsaW5ld2lkdGg9MSxjb2xvdXIgPSAiYmx1ZSIpICsgbGFicyh4ID0gIkRheSIsIHkgPSAiTnVtYmVyIG9mIGFjY2llbnRzIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB0aXRsZSA9IGxldmVscyhkZl90ZXN0JHBvbGljZV9mb3JjZSlbaV0pCmBgYAoKTG9uZG9uOgoKYGBge3J9CmkgPC0gMzAKZGZfdGVzdF9maWx0ZXJlZCA8LSBkZl90ZXN0ICU+JSBmaWx0ZXIocG9saWNlX2ZvcmNlID09IGxldmVscyhwb2xpY2VfZm9yY2UpW2ldKQpwcmVkcyA8LSBwcmVkaWN0KGdhbXNbW2ldXSwgZGZfdGVzdF9maWx0ZXJlZCwgdHlwZSA9ICJyZXNwb25zZSIpCmRmX3Rlc3RfZmlsdGVyZWQgIDwtIGNiaW5kKGRmX3Rlc3RfZmlsdGVyZWQscHJlZHMpICU+JSBhcnJhbmdlKGRheV9vZl95ZWFyKQpgYGAKCmBgYHtyfQpkZl9wbG90IDwtIGRhdGEuZnJhbWUobiA9IGRmX3Rlc3RfZmlsdGVyZWQkbnVtYmVyX29mX2NyYXNoZXMscHJlZCA9IGRmX3Rlc3RfZmlsdGVyZWQkcHJlZHMsIAogICAgICAgICAgICAgICAgICAgICAgeW1heCA9IGRmX3Rlc3RfZmlsdGVyZWQkcHJlZHMgKyBrW2ldKlNbaSxdLAogICAgICAgICAgICAgICAgICAgICAgeW1pbiA9IHBtYXgoIGRmX3Rlc3RfZmlsdGVyZWQkcHJlZHMgLSBrW2ldKlNbaSxdLCByZXAoMCwgMzY1KSksdGQgPSAxOjM2NSkKCmRmX3Bsb3QgJT4lIGdncGxvdCgpICsgZ2VvbV9yaWJib24oYWVzKHg9dGQseSA9IHByZWQseW1pbiA9IHltaW4sIHltYXggPSB5bWF4KSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZmlsbCA9ICJibHVlIiwgYWxwaGEgPSAwLjI1KSArIAogIGdlb21fbGluZShhZXMoeD10ZCx5ID0gbiksbGluZXdpZHRoPTEsY29sb3VyID0gImJsYWNrIikgKyAKICBnZW9tX2xpbmUoYWVzKHg9dGQseSA9IHByZWQpLGxpbmV3aWR0aD0xLGNvbG91ciA9ICJibHVlIikgKyBsYWJzKHggPSAiRGF5IiwgeSA9ICJOdW1iZXIgb2YgYWNjaWVudHMiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHRpdGxlID0gbGV2ZWxzKGRmX3Rlc3QkcG9saWNlX2ZvcmNlKVtpXSkKYGBgCgp0aGUgbW9kZWwgaXMgc3RpbGwgb3QgZ3JhZXQKCg==